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1.  INTRODUCTION 


This  report  describes  the  auroral  spacecraft  charging  model  implemented  in  the  Nascap-2k 
spacecraft  charging  computer  code.  Drs.  Victoria  A.  Davis  and  Myron  J.  Mandell  developed  and 
implemented  the  flux  and  charging  models,  with  some  assistance  from  Gary  A.  Jongeward  with 
the  auroral  environmental  flux  model.  The  use  of  the  Nascap-2k  auroral  charging  model  is 
described  in  Chapter  2.  The  formulas  used  to  describe  the  surface  fluxes  are  given  in  Chapters  3 
through  7.  Chapter  8  compares  the  results  of  sample  calculations  with  those  computed  using 
other  auroral  spacecraft  charging  codes. 


2.  USING  NASCAP-2K 


Choosing  “Polar”  Orbit  and  “Surface  Charging”  Problem  Type  on  the  Problem  tab  (Figure  1) 
specifies  spacecraft  charging  calculations  using  an  auroral  environment.  The  charging  currents 
can  be  computed  using  either  “Analytic  Currents”  (a  strictly  analytic  description  of  the  currents 
described  in  the  sections  Incident  electron  flux  and  Analytic  model  of  incident  ion  flux  below)  or 
“Analytic  Electron  &  Tracked  Ion  Currents”  (where  the  electron  current  is  computed  analytically 
and  the  ion  current  is  computed  by  particle  tracking).  Ion  currents  from  particle  tracking  require 
“Potentials  in  Space”  and  “Surface  Currents”  calculations.  Space  potentials  can  be  computed 
with  either  “Analytic  Space  Charge”  or  “Self-consistent  with  Ion  Trajectories.”  Only  the  low 
energy  component  of  the  plasma  is  used  in  the  computation  of  space  charge,  as  the  high  energy 
electrons  contribute  only  a  low  density  of  order  106  m'3. 
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Figure  1.  Problem  Tab  for  an  Auroral  Charging  Calculation  with  the  Charge  Density  Computed  From 
Tracked  Ions  and  Barometric  Electrons. 


The  components  of  the  environment  are  specified  on  the  Environment  Tab.  (Figure  2) 
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Figure  2.  Environment  Tab  for  Auroral  Charging. 

Only  species  with  positive  charge  are  used  in  the  computation  of  the  tracked  ion  flux  component. 
If  the  space  potentials  are  computed  with  analytic  space  charge,  the  ion  flux  is  computed  using 
the  parameters  given  on  the  “Surface  Currents”  subtab  of  the  “Particles”  tab.  If  the  space 
potentials  are  computed  self-consistently  with  ion  trajectories,  the  ion  flux  is  computed  using  the 
parameters  given  on  the  “Ion  Densities”  subtab  (Figure  3)  of  the  “Particles”  tab. 


:*&MP*m**  :  mmmmA  ******  l 

km Dwitwei'' ■;  / ;  ■ . ■  .  >  •  • ,  •’ 

toHtu*adtoM*xamstat»nOyc<>wii*«lani^  j 


.J 


Figure  3.  Ion  Densities  Subtab  of  the  Particles  Tab. 
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Figure  4.  Advanced  Screen  of  the  Particles  Tab.  Calculations  Self-Consistent  with  Ion  Densities  Require 
Extra  Particles  Created  Either  as  a  Thermal  Distribution  or  by  Subdivision. 

The  initial  potentials  are  specified  in  the  same  way  as  for  geosynchronous  charging  calculations. 


Figure  5.  Initial  Potentials  Tab. 


The  charging  timesteps  are  specified  on  the  charging  tab  (Figure  6)  in  the  same  way  as  for  a 
geosynchronous  charging  calculation.  Computing  surface  charging  self-consistently  with  tracked 
ions  is  an  iterative  calculation.  In  this  case,  the  number  of  iterations  is  given  in  the  Iteration  box 
on  the  Space  Potentials  Tab  (Figure  7)  and  the  number  of  timesteps  on  the  Charging  tab  is 
ignored.  One  timestep  is  used  per  iteration. 
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Figure  6.  Charging  Tab. 
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Figure  7.  Space  Potentials  Tab.  Auroral  Charging  Calculations  Using  Tracked  Ions  Use  the  Iteration 
Parameters  in  the  Space  Potentials  and  Ion  Density(Self-Consistent  Ions  Only)  Calculations. 

The  script  used  to  specify  an  auroral  charging  calculation  with  “Analytic  Currents”  is  the  same 
as  the  script  used  to  specify  a  geosynchronous  charging  calculation  except  for  the  environment 
definition.  The  script  used  to  specify  an  auroral  charging  calculation  with  “Analytic  Electron  & 
Tracked  Ion  Currents”  is  similar  to  that  for  a  LEO  surface  charging  calculation  using  ‘Tracked 
Particle  Currents.”  The  calculation  is  initialized  with  “Charge  Surfaces,”  “Potentials  in  Space,” 
“Create  Particles,”  and  “Track  Particles,”  commands.  The  sequence  is  repeated  “Number  of 
Iterations”  times,  ending  with  a  “Potentials  in  Space”  command.  For  “Self-consistent  with  Ion 
Trajectories”  charge  densities,  the  “Traj_Ion”  command  is  included  in  the  sequence  before  each 


4 


“Potentials  in  Space”  command  after  the  first  one.  The  first  Charge  Surfaces  command  has  the 
same  subcommands  as  an  “Analytic  Currents”  calculation,  except  the  “DoOneTimeStep” 
subcommand  is  used  in  place  of  the  “DoTimeSteps”  subcommand.  A  “UseTrackedlons” 
subcommand  specifies  that  the  most  recent  tracked  surface  current  values  are  to  be  used  as  the 
ion  flux.  A  “FieldsFromFile”  subcommand  specifies  that  the  electric  fields  used  to  compute  the 
suppression  of  secondary  electrons  and  photoelectrons  are  to  be  read  from  the  database  rather 
than  computed  by  the  BEM  module.  The  BEM  module  assumes  no  space  charge,  while  the 
Potentials  in  Space  module  of  Nascap-2k  (which  writes  its  results  into  the  database)  computes 
electric  fields  using  the  specified  model  for  space  charge. 


3.  INCIDENT  ELECTRON  FLUX 


All  species  of  the  low  energy  plasma  are  described  by  a  Maxwellian  distribution.  The  high 
energy  electrons  are  described  by  a  Fontheim  distribution,  which  has  three  components,  specified 
by  the  net  current  in  each  component.  The  Maxwellian  component  describes  a  broad  electron 
distribution,  the  Gaussian  component  describes  the  “inverted- V”  part  of  the  spectrum,  and  the 
Power  Law  component  describes  the  secondary  and  backscattered  electrons  from  interactions 
between  the  Gaussian  beam  and  the  rest  of  the  plasma.  The  Power  Law  component  only 
contributes  at  energies  between  specified  lower  and  upper  cutoffs.  The  effective  density  of  each 
component  is  the  expectation  value  of  the  distribution  function. 

The  electron  differential  flux  (m-2  s-1  eV-1)  is  specified  by 


Flux  (E)  = 


27t0me  0 


nexp 


f 


e 


+  TtCmaxEeXP 


V  emax/ 


+  rcCgauSsEeXP 


f  'E  -E^ 

gauss 


+  < 


power 


E-a 


(1) 


where  n  and  0  are  the  density  and  temperature  of  the  low  energy  ionospheric  plasma,  e  and  me 
are  the  electron  charge  and  mass,  and  the  £s,  0max,  Egauss,  A,  and  a  are  constants.  Figure  8  shows 
the  electron  flux  for  an  auroral  environment. 


Figure  8.  Fontheim  Electron  Spectrum  Using  Parameters  Shown  in  Figure  2. 
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4.  ANALYTIC  MODEL  OF  INCIDENT  ION  FLUX 


The  ion  flux  density  to  each  surface  is  assumed  to  be  the  orbit-limited  ion  flux  density  to  a  point 
on  a  uniform  potential  sphere  at  the  same  angle  relative  to  the  ram  direction.  The  orbit-limited 
assumption,  which  is  appropriate  to  tenuous  plasma  conditions  where  the  Debye  length  is  tens  to 
hundreds  of  meters,  is  questionable  in  the  auroral  case  for  which  the  Debye  length  may  be  of 
order  10  cm.  Also,  the  spherical  geometry  approximation  may  be  less  valid  here  than  for  tenuous 
plasma  environments.  Nonetheless,  use  of  this  model  provides  a  means  of  obtaining  at  least 
qualitative  results  without  the  effort  and  uncertainties  associated  with  particle  tracking. 

We  next  derive  the  equations  for  ion  collection  by  a  negatively  charged  conducting  sphere  of 
radius  a  in  a  flowing  plasma  with  ion  temperature  T.  We  determine  the  current  density  as  a 
function  of  angle.  All  velocities  are  in  units  of  the  plasma  thermal  velocity 

VT  =(kT/M)1/2  (2) 

where  k  is  Boltzman’s  constant  and  M  is  the  ion  mass.  Energies  are  measured  in  units  of  kT,  so 
that  energy  conservation  for  an  ion  at  the  surface  of  an  attracting  sphere  is  expressed  by 

—  v2-<J)  =  E  =  —  Vq  (3) 

2  2 


where  v  is  the  velocity  of  the  ion,  the  potential  energy  of  the  ion  at  the  surface  of  the  sphere,  E 
the  total  energy  of  the  ion,  and  v0  the  speed  of  the  ion  at  an  infinite  distance  from  the  sphere. 

The  normal  current  density  at  a  point  on  the  sphere  where  the  unit  outward  normal  is  n  and  is 
given  by 


j  =  j(n)  = 


v-n  f0(v„)d3v 


(4) 


where  f0  is  the  plasma  velocity  space  distribution  function  and  the  integration  variable  is  the 
negative  of  the  impact  velocity  at  the  surface.  In  the  rest  frame  of  the  sphere 


fo(v0)  =  (2Jt)-3,2exp 


(5) 


where  Vo  is  the  flow  velocity  of  the  plasma. 

Now  introduce  a  coordinate  system  (Figure  9)  whose  z  axis  lies  in  the  direction  of  the  normal  to 
the  sphere  at  the  point  where  we  wish  to  calculate  j.  Let  the  x-z  plane  be  determined  by  this 
normal  and  the  plasma  flow  velocity  and  let  cp  be  the  azimuthal  angle  between  the  x-z  plane  and 
the  orbital  plane  of  an  ion.  The  angle  between  the  surface  normal  and  the  flow  velocity  is  £.  A 
particle  launched  at  the  sphere  surface  with  polar  angle  a  has  polar  angle  0«,  at  r  = 
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Figure  9.  Coordinates  Used  in  the  Drifting  Maxwellian  Ion  Current  Calculation. 


The  current  j  may  now  be  written 

°°  n/2  2n 

j  =  (2tc)_3/2  dvn  |  da 


dtp  v0(vo+20)xexp  ~(vq -2v0V0cos\|/  + V2)Jxcosasina  (6) 


where  a  is  the  polar  angle  of  the  reverse  trajectory  ion  launched  from  the  sphere  and  \)/  is  the 
angle  between  the  flow  vector  and  the  incident  velocity  at  infinity  associated  through  the 
dynamics  with  the  angle  a. 

Referring  to  Figure  9,  we  have 

cos  \|/  =  cos  ^  cos  0„,  +sin£sin0oocos(p  (7) 

Integration  over  the  azimuth  yields 


j  =  (2rc) 


-3/2 


v0(v2  +  2<D)exp(-v2/2)  F(v0)  dv0 


(8) 


where 


n/2 


F(v0)  =  2tt 


exp(-V02 /2)exp(v0V0  cos^cos0£O)xIo  (v0V0  sin  £  sin  0„)  cos  a  sin  a  da  (9) 


and  Io(x)  is  the  modified  Bessel  function  of  zero  order 


1  2n 

I0(x)  =  —  j  exp(-x  cos  cp)d<p 
27t  o 


(10) 
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The  angles  at  infinity  are  functions  of  the  launch  angle  and  the  energy.  The  required  relations  are 
determined  from  the  orbit  relations  in  a  1/r  potential: 


a 

r 


d> 


2(E  +  <I>)sin  a 


-x 


1  +  ^1+— sin  a  cos  (6- 60) 


(ID 


Measuring  orbital  angles  from  the  polar  axis  that  is  normal  to  the  sphere  at  the  launch  point  and 
setting  a  =  r  determines  the  angle  0o, 


cos0o  = 


sin2  a-1 


J 


I  4E(E  +  <t>)  .  2 

1  +  — — ; — -sin  a 
<D2 


— 1 


(12) 


Setting  r  =  <*>  gives 


cos(0^ -0O)  =  -^+—^7— sin2  a  (13) 

The  velocity  and  angle  integrals  are  done  numerically.  Figure  10  shows  the  flux  density  to  a  low- 
potential  sphere,  illustrating  the  progression  from  uniform  flux  a  low  Mach  number  to  cos  0  on 
the  ramward  side  and  zero  on  the  wake  side  for  high  Mach  number.  Figure  1 1  shows  the  flux 
density  to  a  high-potential  sphere,  showing  that  the  flux  decreases  with  ion  energy  (more  or  less 
as  [1-V/E])  and  is  fairly  uniform  over  the  sphere,  even  at  realistic  Mach  numbers. 


Ion  Flux  Density  at  V  =  0.1  Ti 


Figure  10.  Relative  Ion  Currents  at  Low  Sphere  Potentials. 
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Figure  ll.Relative  Ion  Currents  at  About  100  V  in  the  Ionosphere. 


5.  SECONDARY  AND  BACKSCATTERED  ELECTRON  FLUXES 

The  secondary  and  backscattered  electron  fluxes  are  computed  for  each  component  of  the 
incident  current  separately. 

The  secondary  and  backscattered  electron  fluxes  due  to  the  Maxwellian  components  of  the 
incident  electron  fluxes  (low  energy  and  Fontheim)  are  computed  in  the  same  manner  as  for  a 
geosynchronous  Maxwellian  plasma. 

The  secondary  and  backscattered  electron  fluxes  due  to  the  Gaussian  and  Power  Law 
components  of  the  incident  electron  flux  are  computed  by  integrating  the  products  of  the 
distribution  functions  and  the  yield  functions  as  is  done  for  the  Maxwellian  components.  The 
yield  for  the  Maxwellian  and  Power  Law  components  is  taken  to  be  that  for  isotropically  incident 
flux.  The  yield  for  the  Gaussian  component  is  taken  to  be  that  for  normally  incident  flux.  The 
techniques  for  doing  the  integrals  are  given  below. 

The  secondary  electron  flux  due  to  the  incident  ion  flux  when  the  incident  ion  flux  is  computed 
using  the  analytic  orbit  limited  model  is  the  sum  over  the  species  of  the  integral  of  the  product  of 
the  incident  ion  flux  times  the  average  yield  for  a  cosine  theta  distribution  of  ions  with  an  energy 

of  -jMV02  -  e<J) ,  where  M  is  the  ion  mass,  V0  is  the  ram  velocity,  and  (J)  is  the  surface  potential. 

The  rate  of  secondary  emission  is  assumed  to  be  independent  of  mass  (except  indirectly  through 
the  incident  velocity).  If  no  species  are  specified  on  the  Environment  Tab,  the  plasma  is  assumed 
to  be  100  percent  of  ions  of  mass  16  times  the  proton  mass. 
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The  secondary  electron  flux  due  to  the  incident  ion  flux  when  the  incident  ion  flux  is  computed 
by  particle  tracking  is  the  tracked  ion  current  times  the  yield.  The  secondary  yield  is  given  by  a 


normally  incident  ion  with  kinetic  energy 


of  —  MV02  -e<|> ,  where  M  is  16  atomic  mass  units,  V0 
2 


is  the  ram  velocity,  and  <|>  is  the  surface  potential. 


6.  PARTIAL  DENSITIES  AND  COEFFICIENTS 

The  partial  densities  shown  on  the  Environment  Tab  are  given  by  the  expressions  below.  The 
techniques  used  to  compute  the  average  values  of  the  square  root  of  the  inverse  of  the  energy  are 
given  in  the  section  on  integrals  below. 


2  4  |mp  /  / 1 


Density^  =7teCmax0max 


Densitygauss  =  TteC 

gauss  ^  2 


A  exp 


V 

+  E0  Vrt  erfc 


f  E„^ 


V  ^  J 


(E  T  (E2)' 


a-1 


\  power 


(14) 


(15) 


(16) 


The  coefficients  shown  on  the  Environment  Tab  are  given  by  the  expressions  below. 

I 


C  max 


max 

3  max  _ fy.2 


7te0„ 


(17) 


Cgauss  / 

Tie  A—  A  exp 


gauss 


+  E0y/n 


1  +  erf 


(En 


V 


V  A  JJ 


(18) 


=%(«-!) 


7ie 


f  i  1  v 


pa-l  pa-1 


(19) 
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7.  INTEGRALS 


The  technique  used  to  compute  the  average  value  of  the  yield  or  other  quantity  over  the 
distribution  is  given  in  this  section. 

7.1  Average  Value  of  Y(E)  for  Maxwellian 


Y  (E)Eexp 


max  _  max(0,<j>surf ) 


'  idk' 
v  e.  , 


dE 


Eexp 


max(0,<>surf) 


r  ehw'> 


dE 


(20) 


(If  for  ions,  use  <J)surf  rather  than  -(|)surf) 
To  compute  the  numerator  use 

Normalization 


Eexp 


J 


dE 


If  ^surf  >=  0,  then 


=  ej(-01nu  +  <l)surf)du  =  -02}(lnu)du  +  0(t)surf  =02  +  0<|)surf  =  02 


1  + 


fosurf 

0 


If  <l>surf  <  0.  then 


=  -0  exp 


,2 _ j  ^surf 


0 


/(lnu)du  =  02 


(21) 


(22) 


7.2  Average  Value  of  Y(E)  for  Gaussian  Term 


Y(E)  Eexp 


f  ^  ^ gauss  fosurf 


x2\ 


dE 


gauss  __  max(O,0surr) 


{yyauss  _ 


Eexp 


dE 


max(°.«t»surf) 


(23) 


11 


f  {e+e  . (  ,  .  e-e  -<i>surf 

Y(E)Eexp - 8  —  dE  ,  substitute  u  - - - - 


max  (0. ^urf) 


~Epauss  ^surf  *~Egam!i 

A  '  A 


Y  ( Au  +  Egauss  +  (>iurf )  (Au  +  Egauss  +  <t>surf )  exp  (-u2 )  Adu 


Split  into  three  integrals  to  isolate  the  peak. 

if  E  +  d)  <0  the  first  term  in  the  expansion  below  is  zero  and  the  second  term  has  a 

11  Egauss  Tsurf  * 

different  lower  limit,  and  may  also  be  zero. 


7. 

~EgUUSS  fosurf  ^patiss 

A  ’  A 


Y(Au+E„„  +*^)(Au+E_  +  4V,)exp(-u!)Adu 


£ 

+  j  Y(Au  +  Egauss  +  ) ( Au  +  Egauss  +  (|)SUrf ) exp (-u2 ) Adu 
+  J  Y  ( Au  +  Egauss  +<t)smf  )(Au+Egauss  +  (t’surf )  exp  (-u2 )  Adu 


| substitute  x  =  exp(-u2 )  and  take  £  to  be  small) 


12 


Y (E)  Eexp  -  E  -~-ga— 


max(0,(J)j;urf ) 


Y(-AV-lnx+Egauss+4)surf)  -A  + 


E  4"  (b  ^ 

gauss  Tsurf 


Egnuss  “^surf  ^yaus.s 


+  Y(Egauss+(t)Surf)A(EgausS+<l)Surf)^erf(^)  +  f-  Y ( A ( In  X  )  +  Egaus$  +  ^surf  )  A  + 


E  -Mb  , 

gauss  Tsurf 


Normalization 


Eexp  - 


E-Ecaus  -<fc 


gauss  Tsurf  Mg 


max(°’<W) 


Egauss  ^gurf  ^^gguss 

A  ’  A 


""Egauss  0SUrf  E, 


( Au  +  Eganss  +  (J)surf  )exp  (-u2 )  Adu 


uexp(-u2)du  +  A  (Egauss  +  <t>surf ) 


gauss  ^surf  Egauss 

A  ’  A 


'Egauss  ^surf  Egauss 

A  ’  A 


exp(-u2)du 


'Egauss  ftsurf  ~‘EgilUgS 
A  ’  A 


u  exp(-u’  )du  +  A  (E_  +  ♦_  )^erfc('max 


A2  .  Egauss  +  <))SUrf  Egauss  1  ,  A  /V  ,  A  \  _ p„ _ _  f  Egauss  ^surf  Egauss 

yexp  -mm M— ^ +A(EgauSs  +  <Lrf  )“^erfc  max  - - - 


Average  Value  of  Y(E)  for  Power  Law: 


max(E,,E2+<|>surf) 

J  Y(E)E(E-^p-'dE 

/  \p°wer  _  max(0,E)+<i>SUff) _ 

\  /  max(E-.,E-.+0surf ) 

J  E(E-^r-'dE 

max(0,E,  +S>surf ) 


13 


(30) 


max(E2,E2+4w )  / 

J  E(E-^r'dE= 

max(0,E,  +4’«uif )  V 


1  ^max(E2-<t>ljrf,E2)' 


J  dz-^surf 

maxf-^  ,E,)"a+1 


i  \  nm(E2-<|>surf,E2)~0 

1  J 

max(-0KUrf  ,Ej)-a 


=  —  J  ( max  ( E2  ~  ^surf  ^  E2  )_“+I  -  max  (-<)£urf ,  E, )  “+1) 

H>,urf  [~](max(E2 -0surf,E2r-ImxKurf,Eir) 


dz 


8.  VALIDATION  OF  RESULTS 

8.1  Comparison  with  SEE  Spacecraft  Charging  Handbook 

The  SEE  Spacecraft  Charging  Handbook  gives  the  following  results  for  a  2  m  Kapton  cube  with 
orbital  speed  of  7500  m/s  in  the  x  direction.  The  charging  time  is  10  s  with  40  timesteps.  The 
environment  is  the  “DMSP1”  environment  shown  in  Figure  13. 


j  4B33P  =  ;  Spacecraft  Name:  |  unnamed  3  F  Eclipse  Environment:  fOMSPi* 


33  Charging  Time  (s):  pfo  Timesteps:  [40 


Solar  array  angle  (deg): 


Orbital  speed  (m/s): 
[7500 


Ram  direction: 

F  |B  |o 


Material: 

Color: 

Dark  Are  a  (m2) 

Max  Potential  (kV) 
Min.  Potential  (kV) 
Max  Differential  (KV) 
Min  Differential  (kV) 
Initial  Current  (A) 
Final  Current  (A) 


Kapton 

|  Red  3 

24  0 

-9  342e-3 
-0.04336 
0  010421 
-0  023602 
-1  22e-3 
6  06e-9 


Figure  12.  Aurora!  Charging  Test  Case  as  Computed  by  SEE  Spacecraft  Charging  Handbook. 
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fitvt*vnm*M4s 


\ 

DMSP1 

ionosphere  Density  (m*3): 

|3.55e9 

Ionosphere  Temperature  (eV): 

|02 

Fraction  Hydrogen: 

Maxwellian  Current  (Am'2): 

|l,40*-6 

Maxwellian  Temperature  (eV): 

[8000  ” 

Maxwellian  coeffieent  (Am'2). 

4  35e4 

Maxwellian  Density  (m'3): 

5.8185 

Gaussian  Current  (Am-2): 

Jl  40e-5 

Gaussian  Energy  (eV) 

[24000 

Gaussian  Width  (eV): 

[16000 

Gaussian  Coeffieent  (Am'2). 

4  0Be4 

Gaussian  Density  (m-3): 

7.77e5 

Power  Law  Current  (Am'2): 

1 670*7 

First  Power  Law  Energy  (eV): 

[50 

Second  Power  Law  Energy  (evy 

Power  Law  Exponent 

In 

Power  Law  Coeffieent  (Am^eV*1): 

3  05et1 

Power  Law  Den  s  Ity  (m*3): 

1  03e6 

Total  Electron  Current  (Am'2): 

-5  B6e-5 

Figure  13.  Environment  Used  for  Test  Cases. 

With  a  minimum  timestep  of  0.01  and  a  maximum  timestep  of  1  s,  the  Nascap-2k,  Polar, 
Analytic  Currents,  charging  model  gives  the  same  results.  For  these  parameters,  this  approach 
gives  too  little  charging.  The  ion  currents  are  overestimated  as  Debye  screening  is  ignored. 


Figure  14.  Results  From  Nascap-2k  Test  Case  Using  Environment  Currents. 

Doing  the  same  problem,  but  tracking  particles  in  from  a  sheath  (located  at  -0.277  V)  gives  more 
charging  as  expected,  since  the  orbit  limited  approximation  overestimates  the  ion  current.  Three 
grids  were  used.  The  main  grid  is  24  cubed  with  a  mesh  size  of  0.8  m.  The  second  grid  is  also  24 
cubed  filling  the  inner  eighth  of  the  main  grid.  The  innermost  grid  is  16  cubed  and  fills  the  center 
of  the  grid.  The  results  are  shown  in  Figure  15.  The  highest  potential  surfaces  are  along  the  front 
edge  of  the  side  and  the  lowest  in  the  back.  Figure  16  shows  the  potential  contours  around  the 
cube.  The  cube  is  moving  to  the  right.  Most  of  the  sheath  is  in  the  second  grid,  which  has  a 
sheath  potential  of  -0.277  V. 
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Figure  15.  Results  From  Nascap-2k  Test  Case  Using  Tracked  Ion  Currents.  The  Top  Left  Surface  is  on  the 
Wake  Side. 


SI 


« 


Figure  16.  Potential  Contours  Around  Charged  Cube  with  Different  Scales. 


Repeating  the  calculation,  but  computing  the  charge  density  using  full  trajectory  ions,  gives 
somewhat  lower  potentials.  Ions  with  a  thermal  distribution  (1.0  in  x;  0.3,  0.4,  0.3  in  y;  and  0.3, 
0.4,  0.3  in  z)  are  tracked  from  the  problem  boundaries.  The  results  are  shown  in  Figure  17.  The 
potentials  are  similar,  although  more  negative.  The  highest  potential  is  at  the  back  edge  of  the 
side. 


Figure  17.  Results  From  Nascap-2k  Test  Case  Using  Tracked  Ion  Currents  and  Self-consistent  Charge 
Density.  The  Top  Left  Surface  is  on  the  Wake  Side. 


There  are  two  significant  differences  between  the  calculation  shown  in  Figure  15  and  that  shown 
in  Figure  17.  The  macroparticles  tracked  are  created  differently  and  the  charge  density  is 
computed  differently.  In  the  first  case,  mono-energetic  ions  are  tracked  from  a  sheath  boundary 
at  -0.277  V,  while  in  the  second  case,  a  thermal  distribution  is  tracked  from  the  problem 
boundaries.  In  the  first  case  the  charge  density  is  computed  using  an  analytic  formulation  and  in 
the  second  the  charge  density  is  computed  from  the  tracked  ions  and  a  barometric  relation  for  the 
electrons.  To  separate  these  two  effects  two  additional  calculations  were  done.  In  the  first,  the 
charge  density  was  computed  using  the  analytic  formulation  and  the  current  computed  from  a 
thermal  distribution  of  ions  tracked  from  the  problem  boundaries.  In  the  second,  the  charge 
density  was  computed  using  ions  tracked  from  the  problem  boundary,  and  the  current  was 
computed  from  ions  tracked  from  a  sheath  at  -0.277  V.  The  potential  time  histories  are  shown  in 
Figure  18.  The  tracked  ion  current  collected  is  shown  in  Figure  19.  In  all  cases,  the  collected 
current  is  roughly  linear  with  chassis  potential.  The  highest  current  gives  the  least  charging.  Both 
the  charge  density  and  the  initial  macroparticle  distribution  influence  the  results.  For  this 
particular  problem,  it  is  not  clear  which  approach  is  more  correct. 
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Figure  18.  Time  History  of  Potentials  for  Different  Charge  Density  Formulations  and  Macroparticle  Initial 
Conditions. 


Figure  19.  Collected  Current  as  a  Function  of  Chassis  Potential  (V)  for  Different  Charge  Density 
Formulations  and  Macroparticle  Initial  Conditions. 


At  lower  densities,  the  analytic  estimation  of  the  ion  current  should  be  more  accurate.  We  did  the 
same  calculation  with  an  ambient  density  two  orders  of  magnitude  lower.  It  is  necessary  to  use  a 
larger  grid,  so  that  the  entire  sheath  fits.  Changing  the  outer  grid  mesh  unit  to  5  meters,  but 
leaving  the  grid  structure  the  same,  gives  a  large  enough  grid.  The  sheath  edge  is  partially  in  the 
outer  grid  and  partially  in  the  second  grid.  The  appropriate  sheath  potentials  are  0.372465  V  and 
0.147813  V.  The  calculations  were  done  using  the  average  sheath  potential  of  0.277  V.  At  the 
lower  density,  the  charging  occurs  faster,  so  we  used  50  timesteps,  with  a  minimum  of  0.001s 
and  a  maximum  of  1  s.  The  final  potentials  are  shown  in  Table  1. 
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Table  1.  Final  Potentials  in  Low  Density  Case. 


Case 

Potential  at  10  s  (V) 

Ions 

Chassis 

Maximum 

Minimum 

None 

-1257 

-1238 

-1313 

Sheath 

■sm 

-997 

■ni 

Boundary 

Self-consistent 

i 

In  the  Sheath/ Analytic  space  charge  case,  the  most  negative  potentials  are  on  the  wake  side  and 
the  least  negative  potentials  are  on  the  ram  side.  In  the  Boundary/Self-consistent  space  charge 
case,  the  most  negative  potentials  are  along  the  side  toward  the  ram  and  the  least  negative 
potential  is  on  the  center  surface  of  the  wake  side. 

Figure  20  shows  sample  trajectories  for  the  Sheath/ Analytic  plasma  case.  Of  the  trajectories 
shown,  the  outer  four  are  undisturbed,  the  inner  two  are  immediately  collected,  and  the  rest  loop 
around  the  cube  one  or  more  times  before  leaving  the  computational  space. 


Figure  20.  Sample  Trajectories  for  Sheath/Analytic  Space  Charge  Case  (Low  Density). 

Figure  21  shows  sample  trajectories  for  the  Boundary/Self-consistent  plasma  case.  Of  the  trajectories 
shown,  the  outer  six  are  undisturbed,  the  inner  four  are  immediately  collected,  the  seventh  from  the 
top  circles  a  couple  of  times  are  then  leaves  the  computational  space,  and  the  eighth  from  the  top 
provides  most  of  the  tangle  of  trajectory  seen  before  being  collected.  The  trajectories  that  orbit  about 
the  cube  before  being  collected  or  leaving  the  computational  space  provide  Debye  shielding  that 
decreases  the  ion  current,  increasing  the  charging.  The  algorithm  for  computing  self-consistent 
potentials  can  give  potentials  that  are  not  as  smoothly  varying  as  desired.  Better  results  are  often 
obtained  by  sharing  the  potentials  or  ion  densities  from  iteration  to  iteration.  The  parameters  for  this 
sharing  are  set  on  the  Space  Potentials  Tab.  In  this  case,  the  space  potentials  have  some  non-physical 
structures  that  may  be  numerically  reducing  the  ion  energies,  leading  to  higher  current  collection  than 
physically  correct. 
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Figure  21.  Sample  Trajectories  for  Boundary/Seif-Consistent  Case  (Low  Density). 


8.2  Comparison  with  POLAR  Results 

To  validate  Nascap-2k  we  compare  the  results  obtained  using  Nascap-2k  with  those  obtained 
using  POLAR,  which  has  been  validated.  The  POLAR  calculations  to  which  we  will  compare  are 
given  in  Reference  1. 

We  used  the  POLAR  object,  which  generates  some  “inconsistent  edge”  error  messages,  but  is 
valid.  To  create  the  object,  we  used  Nascap/GEO  for  Windows  to  create  the  object,  read  the 
Nascap/GEO  object  into  DynaPAC,  and  then  read  the  DynaPAC  object  into  Object  Toolkit. 


Figure  22.  POLAR  DMSP  Object. 


Charging  in  eclipse  from  an  initial  potential  of  -10  V  was  computed  for  a  period  of  6  s  with  45 
timesteps  varying  from  0.0125  s  to  0.2  s.  In  all  cases,  the  charging  currents  were  computed  using 
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an  analytic  model  for  the  electrons  and  tracked  ions.  The  charge  density  was  computed  either 
analytically  or  self-consistently  with  the  ion  trajectories.  In  the  first  case,  the  ion  currents  are 
computed  by  tracking  macroparticles  from  the  sheath  and  in  the  second  case  by  tracking 
macroparticles  from  the  boundary  of  the  computational  space.  The  environment  used  is  shown  in 
Figure  23. 
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Figure  23.  Environment  Used  for  DMSP  Charging  Calculations. 


The  results  of  the  calculations  are  given  in  Figure  24  through  Figure  31.  With  these  choices  for 
parameters,  the  resulting  chassis  potential  for  the  analytic  charge  density  case  is  -707  V  and  for 
the  self-consistent  charge  density  case  is  -845  V.  In  the  analytic  charge  density  case  the 
potentials  vary  between  -414  and  -888  V;  while  in  the  self-consistent  charge  density  case  the 
potentials  vary  between  -710  and  -1025  V.  The  distribution  of  surface  potentials  is  essentially  the 
same  for  the  two  cases,  with  well-shadowed  surfaces  the  most  negative.  The  chassis  rapidly 
charges  to  the  -400  to  500  V  level  and  continues  to  charge  at  about  50  to  60  V  per  second.  With 
more  timesteps,  the  results  might  be  slightly  different. 
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Figure  24.  Resulting  Surface  Potentials  for  Case  with  Analytic  Charge  Density. 


Figure  25.  Resulting  Surface  Potentials  for  Case  with  Self-Consistent  Charge  Density. 
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Analytic  charge  density 


Figure  26.  Time  History  of  Charging  with  Analytic  Charge  Density. 


Self-consistent  charge  density 


Figure  27.  Time  History  of  Charging  with  Self-Consistent  Charge  Density. 


Several  interesting  features  can  be  seen  in  the  space  potentials.  In  the  analytic  space  charge  case, 
the  potentials  follow  the  contours  of  the  satellite  more  closely  and  the  influence  of  the  wake  in 
the  higher  potential  regions  is  more  pronounced.  In  the  self-consistent  space  charge  case,  the 
potential  falls  off  a  bit  faster  in  the  wake,  and  the  sheath  edge  is  much  smoother  on  the  ram  side. 
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Figure  28.  Resulting  Space  Potentials  for  Case  with  Analytic  Charge  Density. 


Figure  29.  Resulting  Space  Potentials  for  Case  with  Self-Consistent  Charge  Density. 
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Figure  30.  Resulting  Sheath  Location  for  Case  with  Analytic  Charge  Density. 


Figure  31.  Resulting  Sheath  Location  for  Case  with  Self-Consistent  Charge  Density. 


These  results  can  be  compared  with  Figure  5  of  Reference  1,  shown  below  as  Figure  32.  The 
Nascap-2k  results  do  not  have  the  presumably  spurious  hump  at  0.2  sec.  The  Nascap-2k  results 
show  more  charging  and  are  continuing  to  charge.  The  relative  charging  of  the  different  surfaces 
are  somewhat  different.  As  the  selected  surfaces  are  probably  different,  no  firm  conclusions  can 
be  drawn. 
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Seconds 


Figure  32.  Figure  5  of  Reference  1. 

We  also  repeated  another  set  calculations  in  the  same  paper.  We  set  the  thickness  of  Teflon  to 
2.8e-3  and  changed  all  the  Kapton  to  Teflon.  The  results  are  shown  in  Figure  33  through  Figure 
35.  These  results  can  be  compared  with  Figure  36  through  Figure  38  obtained  from  POLAR. 
Again,  the  results  using  the  two  codes  are  similar  in  character.  The  shaded  surfaces  are  the  most 
negative  and  the  ram-wake  difference  is  small.  Using  POLAR,  the  wake  side  charges  more  than 
the  ram  side.  Using  Nascap-2k,  the  ram  side  charges  slightly  more  than  the  wake  side.  The 
incident  ions  are  focused  onto  the  wake  side  of  the  spacecraft. 
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Figure  33.  Resulting  Surface  Potentials  for  Case  with  All  Teflon  and  Self-Consistent  Charge  Density. 
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Figure  34.  Resulting  Surface  Potentials  for  Case  with  All  Teflon  and  Self-Consistent  Charge  Density 


Self-consistent  charge  density,  other  example 


Time  (s) 

Figure  35.  Time  History  of  Charging  with  All  Teflon  and  Self-Consistent  Charge  Density. 
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Figure  36.  Resulting  Surface  Potentials  Using  POLAR  for  Teflon  Only  Case  From  Reference  1. 
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Figure  37.  Resulting  Surface  Potentials  Using  POLAR  for  Teflon  Only  Case  From  Reference  1. 


Figure  38.  Time  History  of  Surface  Potentials  Using  POLAR  for  Teflon  Only  Case  From  Reference  1. 
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